Laplacian smoothing

library(plotly)

## Warning: package 'plotly' was built under R version 3.3.3

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(microbenchmark)

## Warning: package 'microbenchmark' was built under R version 3.3.3

grid <- read.csv("~/GitHub/SDS385-course-work/Excercise 8 Spatial smoothing/Laplacian smoothing/fmri_z.csv")
grid=as.matrix(grid)
source('~/GitHub/SDS385-course-work/Excercise 8 Spatial smoothing/Laplacian smoothing/functions.R')
p1=plot_ly(z=grid,type='heatmap')

n=nrow(grid)
y=as(c(grid),"CsparseMatrix")

lambda=0.1
C=makeC_sparse(n,n,lambda)
x = solve(C,y,sparse=TRUE)
smoothed_grid = matrix(data=x,nrow=n)
p2=plot_ly(z=smoothed_grid,type='heatmap')

lambda=1
C=makeC_sparse(n,n,lambda)
x = solve(C,y,sparse=TRUE)
smoothed_grid = matrix(data=x,nrow=n)
p3=plot_ly(z=smoothed_grid,type='heatmap')

lambda=10
C=makeC_sparse(n,n,lambda)
x = solve(C,y,sparse=TRUE)
smoothed_grid = matrix(data=x,nrow=n)
p4=plot_ly(z=smoothed_grid,type='heatmap')

p=subplot(p1,p2,p3,p4,nrows=2)
p

lambda=1
C=makeC_sparse(n,n,lambda)

result=gauss_seidel(C,y,50,targetfunction = TRUE)
targetfunction=result[[2]]
plot(targetfunction,type='l',col='blue')
legend("topright",col=c('blue','red','black'),legend=c("Gauss Seidel","Jacobi","Direct solver"),lty=1)
result=jacobi(C,y,50,targetfunction = TRUE)
targetfunction=result[[2]]
lines(targetfunction,col='red')
x = solve(C,y,sparse=TRUE)
direct_solver=target_function(C,y,x)[1,1]
lines(rep(direct_solver,50),col='black')

x=gauss_seidel(C,y,8)[[1]]
smoothed = matrix(data=x,nrow=n)
plot_ly(z=smoothed,type='heatmap')

x=jacobi(C,y,15)[[1]]
smoothed = matrix(data=x,nrow=n)
plot_ly(z=smoothed,type='heatmap')

microbenchmark(gauss_seidel(C,y,8),times=1)$time

## [1] 28099008

microbenchmark(jacobi(C,y,15),times=1)$time

## [1] 31740398

microbenchmark(solve(C,y,sparse=TRUE),times=1)$time

## [1] 4680716